#pragma once

#include "Vector.hpp"

template <class T>
class Matrix : public Vector<T>
{
protected:
    const int Row, Col;

public:
    Matrix(const int Row, const int Col, const T &t = T()) : Row(Row), Col(Col), Vector<T>(Row * Col, t) {} // 初始化构造函数
    Matrix(const int Row, const int Col, std::initializer_list<T> l) : Row(Row), Col(Col), Vector<T>(l) {}  // 通过列表初始化
    Matrix(const int Row, const int Col, const T arr[]) : Row(Row), Col(Col), Vector<T>(Row * Col, arr) {}  // 通过列表初始化

    // 一定要注意还有普通拷贝构造
    Matrix(const Matrix &rhs) : Row(rhs.row()), Col(rhs.col()), Vector<T>(rhs.size())
    {
        for (int d = 0; d < this->Dim; d++)
            this->value[d] = rhs[d];
    }

    // 只能显式地调用此构造函数
    template <class T2>
    explicit Matrix(const Matrix<T2> &rhs) : Row(rhs.row()), Col(rhs.col()), Vector<T>(rhs.size())
    {
        for (int d = 0; d < this->Dim; d++)
            this->value[d] = static_cast<T>(rhs[d]);
    }

    // 返回单位阵
    static Matrix<T> unit(int row)
    {
        Matrix<T> I(row, row);
        for (int d = 0; d < row; d++)
            I(d, d) = 1;
        return I;
    }

    // 注意在模板类中无法直接访问基类成员，需要通过 this 指针调用
    T &operator()(int i, int j) { return this->value[i * Row + j]; }
    T operator()(int i, int j) const { return this->value[i * Row + j]; }

    // 返回第 j 列
    Vector<T> operator()(int j) const
    {
        Vector<T> v(Row);
        for (int i = 0; i < Row; i++)
        {
            v[i] = this->value[i * Row + j];
        }
        return v;
    }

    // 在 r,c 位置插入列向量 v
    template <class T2>
    void insert(int r, int c, Vector<T2> &v)
    {
        int maxRow = MIN(v.size(), Row - r);
        for (int i = r; i < maxRow + r; i++)
        {
            this->value[r * Row + c] = v[i - r];
        }
    }

    // 获取矩阵的子阵
    Matrix<T> subMatirx(Vector<int> &horizon, Vector<int> &vertical)
    {
        int rowNum = horizon.size(), colNum = vertical.size();
        Matrix<T> sub(rowNum, colNum);
        for (int i = 0; i < rowNum; i++)
        {
            for (int j = 0; j < colNum; j++)
            {
                sub(i, j) = this->value[horizon[i] * Row + vertical[j]];
            }
        }
        return sub;
    }

    int row() const { return Row; }
    int col() const { return Col; }
    int order() const { return MIN(Row, Col); }
    bool square() const { return (Row == Col ? true : false); }

    friend std::ostream &operator<<(std::ostream &out, const Matrix<T> &mat)
    {
        int Row = mat.row(), Col = mat.col();
        for (int i = 0; i < Row; i++)
        {
            for (int j = 0; j < Col; j++)
            {
                out << std::setw(12) << mat(i, j) << " ";
            }
            out << std::endl;
        }
        return out;
    }

// 加一个矩阵/在对角线上加一个向量
#define MATRIX_BINARY_OP(OpName, Op)                       \
    template <class T2>                                    \
    auto OpName(const Matrix<T2> &rhs) const               \
    {                                                      \
        using Tx = decltype(this->value[0] Op rhs[0]);     \
        Matrix<Tx> res(Row, Col);                          \
        for (int i = 0; i < this->Dim; i++)                \
            res[i] = this->value[i] Op rhs[i];             \
        return res;                                        \
    }                                                      \
    template <class T2>                                    \
    auto OpName(const Vector<T2> &rhs) const               \
    {                                                      \
        using Tx = decltype(this->value[0] Op rhs[0]);     \
        Matrix<Tx> res(Row, Col, this->value);             \
        for (int i = 0; i < rhs.size(); i++)               \
            res[i * Row + i] = res[i * Row + i] Op rhs[i]; \
        return res;                                        \
    }

    MATRIX_BINARY_OP(operator+, +)
    MATRIX_BINARY_OP(operator-, -)
#undef MATRIX_BINARY_OP

// 矩阵赋值（定义矩阵以后，可以直接用向量赋值）
#define MATRIX_EVALUATE(Type)               \
    void operator=(const Type &rhs)         \
    {                                       \
        for (int i = 0; i < this->Dim; i++) \
        {                                   \
            this->value[i] = rhs[i];        \
        }                                   \
    }

    MATRIX_EVALUATE(Vector<T>)
    MATRIX_EVALUATE(Matrix)
#undef MATRIX_EVALUATE

    // 乘除法
#define MATRIX_BINARY_OP(OpName, Op)        \
    auto OpName(const T &rhs) const         \
    {                                       \
        Matrix res(Row, Col);               \
        for (int i = 0; i < this->Dim; i++) \
            res[i] = this->value[i] Op rhs; \
        return res;                         \
    }

    MATRIX_BINARY_OP(operator*, *)
    MATRIX_BINARY_OP(operator/, /)
#undef MATRIX_BINARY_OP

    // 矩阵乘法
    Vector<T> operator*(const Vector<T> &_v)
    {
        Vector<T> v(Row);
        T sum;
        for (int i = 0; i < Row; i++)
        {
            sum = 0;
            for (int j = 0; j < Col; j++)
            {
                sum += operator()(i, j) * _v[j];
            }
            v[i] = sum;
        }
        return v;
    }
    Matrix operator*(const Matrix &_m)
    {
        int _column = _m.col();
        Matrix m(Row, _column);
        T sum;

        for (int i = 0; i < Row; i++)
        {
            for (int j = 0; j < _column; j++)
            {
                sum = 0;
                for (int k = 0; k < Col; k++)
                {
                    sum += operator()(i, k) * _m(k, j);
                }
                m(i, j) = sum;
            }
        }
        return m;
    }
};

// ====================================================================================================

struct Mat
{
    template <class T>
    static bool approx(const Matrix<T> &lhs, const Matrix<T> &rhs)
    {
        if (rhs.row() != lhs.row() || rhs.col() != lhs.col())
        {
            return false;
        }
        for (int i = 0; i < lhs.row(); i++)
        {
            for (int j = 0; j < lhs.col(); j++)
            {
                if ((lhs(i, j) - rhs(i, j)) > ACCURACY)
                {
                    return false;
                }
            }
        }
        return true;
    }

    // 置换作用于矩阵
    template <class T>
    static Matrix<T> exchange(const Perm &perm, const Matrix<T> &A)
    {
        int row = A.row(), col = A.col();
        Matrix<T> pA(row, col);
        for (int i = 0; i < row; i++)
        {
            for (int j = 0; j < col; j++)
            {
                pA(i, j) = A(perm[i], j);
            }
        }
        return pA;
    }

    static Matrix<Real> Hilbert(int n)
    {
        Matrix<Real> A(n, n);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                A(i, j) = 1.0 / (i + j + 1);
            }
        }
        return A;
    }

    // 多项式的拉普拉斯阵
    static Matrix<Real> Laplace(int n, Real a[])
    {
        Matrix<Real> A(n, n);
        for (int i = 0; i < n; i++)
        {
            A(0, i) = a[i];
            if (i < n - 1)
            {
                A(i + 1, i) = -1;
            }
        }
        return A;
    }

    static Matrix<Real> rand_mat(int row, int column)
    {
        Matrix<Real> m(row, column);
        srand((int)time(0));
        for (int i = 0; i < row; i++)
        {
            for (int j = 0; j < column; j++)
            {
                m(i, j) = rand() / Real(RAND_MAX);
            }
        }
        return m;
    }

    // ====================================================================================

    template <class T>
    static T norm(const Matrix<T> &lhs, int nt = 2)
    {
        int Row = lhs.row(), Col = lhs.col();
        if (nt == 2)
        {
            return sqrt(Vec::dot(lhs, lhs));
        }
        else if (nt == 1)
        {
            T max = 0;
            for (int i = 0; i < Col; i++)
            {
                T a = 0;
                for (int j = 0; j < Row; a += fabs(lhs(j++, i)))
                    ;
                max = (a > max ? a : max);
            }
            return max;
        }
        else if (nt == 0)
        {
            T max = 0;
            for (int i = 0; i < Row; i++)
            {
                T a = 0;
                for (int j = 0; j < Col; a += fabs(lhs(i, j++)))
                    ;
                max = (a > max ? a : max);
            }
            return max;
        }
        return 0; // never reached
    }

    // e 范数 / f 范数
    template <class T>
    static T norm(const Matrix<T> &lhs, char c)
    {
        int Row = lhs.row(), Col = lhs.col();
        if (c == 'e')
        {
            T sum = 0;
            for (int i = 0; i < Row; i++)
            {
                for (int j = 0; j < Col; j++)
                {
                    if (i != j)
                    {
                        sum += lhs(i, j) * lhs(i, j);
                    }
                }
            }
            return sqrt(sum);
        }
        else if (c == 'f')
        {
            T sum = 0;
            for (int i = 0; i < Row; i++)
            {
                for (int j = 0; j < Col; j++)
                {
                    sum += lhs(i, j) * lhs(i, j);
                }
            }
            return sqrt(sum);
        }
        return 0;
    }

    // get upper or lower
    template <class T>
    static Matrix<T> upper(Matrix<T> &A, bool hasDialog = true)
    {
        int Col = A.col(), Row = A.col();
        Matrix<T> m(Row, Col);
        int limit = MIN(Row, Col);
        for (int i = 0; i < limit; i++)
        {
            for (int j = (hasDialog ? i : i + 1); j < limit; j++)
            {
                m(i, j) = A(i, j);
            }
        }
        return m;
    }

    template <class T>
    static Matrix<T> lower(Matrix<T> &A, bool hasDialog = true)
    {
        int Col = A.col(), Row = A.col();
        Matrix<T> m(Row, Col);
        int limit = MIN(Row, Col);
        for (int i = 0; i < limit; i++)
        {
            for (int j = 0; j < (hasDialog ? i + 1 : i); j++)
            {
                m(i, j) = A(i, j);
            }
        }
        return m;
    }

    template <class T>
    static Vector<T> diag(Matrix<T> &A)
    {
        int Col = A.col(), Row = A.col();
        int limit = MIN(Row, Col);
        Vector<T> v(limit);
        for (int i = 0; i < limit; i++)
        {
            v[i] = A(i, i);
        }
        return v;
    }

    template <class T>
    static Matrix<T> transpose(Matrix<T> &A)
    {
        int Col = A.col(), Row = A.col();
        Matrix<T> m(Col, Row);
        for (int i = 0; i < Col; i++)
        {
            for (int j = 0; j < Row; j++)
            {
                m(i, j) = A(j, i);
            }
        }
        return m;
    }

    template <class T>
    static Matrix<T> power(Matrix<T> &A, int p)
    {
        Matrix<T> P = Matrix<T>::unit(A.row());
        for (int i = 0; i < p; i++)
        {
            P = P * A;
        }
        return P;
    }

    template <class T>
    static T det(Matrix<T> &A)
    {
        int sgn = -1, n = A.row();
        if (n == 1)
        {
            return A[0];
        }
        T sum = 0;
        Vector<int> rowIndex(n - 1), colIndex(n - 1);
        // 子矩阵的行列向量标记
        for (int i = 1; i < n; i++)
        {
            rowIndex[i - 1] = i;
            colIndex[i - 1] = i;
        }
        Matrix<T> B(n - 1, n - 1);
        for (int i = 0; i < n; i++)
        {
            sgn *= -1;
            B = A.subMatirx(rowIndex, colIndex);
            sum += sgn * A[i] * det(B);
            colIndex[i] = i;
        }
        return sum;
    }

    // template <class T>
    // Complex eigen2(const Matrix<T> &A, int row = 0, int col = 0)
    // {
    //     T b = -1 * (A(row, col) + A(row + 1, col + 1));
    //     T c = A(row, col) * A(row + 1, col + 1) - A(row, col + 1) * A(row + 1, col);
    //     Complex z = sqrt_C(b * b - 4 * c);
    //     z = (z - b) / 2;
    //     return z;
    // }
};